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The problem of construction of the Wannier functions (WFs) in a restricted Hilbert space of 
eigenstates of the one-electron Hamiltonian H (forming the so-called low-energy part of the spec- 
trum) can be formulated in several different ways. One possibility is to use the projector-operator 
techniques, which pick up a set of trial atomic orbitals and project them onto the given Hilbert 
space. Another possibility is to employ the downfolding method, which eliminates the high-energy 
part of the spectrum and incorporates all related to it properties into the energy-dependence of 
an effective Hamiltonian. We show that by modifying the high-energy part of the spectrum of the 
original Hamiltonian H, which is rather irrelevant to the construction of WFs in the low-energy part 
of the spectrum, these two methods can be formulated in an absolutely exact and identical form, so 
that the main difference between them is reduced to the choice of the trial orbitals. Concerning the 
latter part of the problem, we argue that an optimal choice for trial orbitals can be based on the 
maximization of the site-diagonal part of the density matrix. This idea is illustrated for a simple toy 
model, consisting of only two bands, as well as for a more realistic example of t^g bands in V2O3. 
Using the model analysis, we explicitly show that a bad choice of trial orbitals can be linked to the 
discontinuity of phase of the Bloch waves in the reciprocal space, which leads to the derealization 
of WFs in the real space. Nevertheless, such a discontinuity does not necessary contribute to the 
matrix elements of H in the Wannier basis. Similar tendencies are seen in realistic calculations for 
V2O3, though with some variations caused by the multi-orbital effects. An analogy with the search 
of the ground state of a many-electron system is also discussed. 

PACS numbers: 71.15.Ap, 71.10.Fd, 71.28.+d, 71. 15. Mb 



I. INTRODUCTION 



The Wannier functions (WFs) play a key role in the solid-state physics, as they make a direct connection between 
the reciprocal k-momentum space and the real space representations for the quantum- mechanical operators, which 
allows us to formulate the problem of electronic-structure calculations for the infinite periodical systems in a basis 
of localized (Wannier) orbitalsii The construction of WFs using first-principles electronic structure calculations has 
attracted a considerable attention re centlyi2iiMi£i£iL2i2ii2iiI The reason is mainly twofold: 

(a) WFs appear to be extremely useful for the visualization and interpretation of the behavior of the local quantities 
in the real space. For example, in the electronic structure calculations based on the plane waves, such a 
construction can add many new functionalities, which were initially a privilege of the methods working in the 
basis of localized atomic orbitalsiSiiS 

(b) Another important application of WFs is related with the construction of the so-called " ab initio models for 
the strongly-correlated materials" , which are designed to deal with the low-energy properties of these systems 
and go beyond the conventional (and oversimplified in the case of the strongly-correlated systems) local-density 
approximation (LDA). In this case, the Wannier orbitals form the basis of the physical Hilbert space, which is 
then used for the construction of the model (typically, Hubbard-type) Hamiltoniansii^i£*iiiL2ii2iii 

The choice of the Wannier function is not unique, and their extension in the real space can vary depending on the 
procedure and the parameters of calculations. In order to illustrate the importance of this problem, one may remind 
to the reader that one extreme example of WFs are the Bloch states, which are totally delocalized in the real space. 
Of course, the physics should not depend on whether the WFs are localized or not, because all representations provide 
the complete sets of the basis functions (for the low-energy part of the spectrum), which are totally equivalent from 
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the physical point of view. However, from the viewpoint of practical applications, it is more convenient to work with 
the localized orbitals. For example, in the process of construction of the Hubbard Hamiltonian one would always like 
to keep only the site-diagonal elements of Coulomb and exchange interactions and neglect all intersite interactions. 
This is only possible when the Wannier basis is sufficiently well localized in the real space, so that the overlap between 
the functions centered at different lattice sites is minimal. The localization of the basis functions is also important in 
applications dealing with the response of the electronic system onto a perturbation confined within a single atomic 
site, like the orbital magnetization in solids^ 

Since the idea of construction of "ab initio models for the strongly-correlated systems" is still relatively new, there 
is no clear consensus among researchers working in this area about the origin of apparent differences, accuracy, and 
the scopes of applicability of their results i^&SiiSiii We believe that in such a situation it is very important to make 
a direct comparison between different methods of calculation of WFs, and prove several theorems regarding both the 
formulation of the methods and the choice of the computational parameters, which directly control the localization 
of WFs. Particularly, we will prove rigorously that the downfoldingmethod, considered in Refs. 03, El andllll can be 
reformulated as a projector-operator method, employed in Refs.UQ, and0 (and, as an initial stage, in Ref . 0) ■ This 
can be done by a simple (scissor-operator-like) change of the full Hamiltonian of the system in the high-energy part 
of the s pect rum, which does not contribute to the low-energy part. This procedure has been already implemented in 
Refs. IH Tlol and Therefore, two methods are simply equivalent, and the main difference comes from the form of 
trial orbitals, which are used to generate the Wannier function in the projector-operator method and, directly, the 
effective low-energy Hamiltonian in the downfolding method. Concerning the last part of the problem, we will argue 
that an optimal choice for the trial wavefunctions can be based on the maximization of the site-diagonal part of the 
density matrix, as it was originally proposed in Ref. 0- We will also consider an analogy with the search of the ground 
state of the many-electron systems. 



II. ATOMIC ORBITALS AND WANNIER FUNCTIONS 



We assume that there is certain set of orthonormalized atomic-like orbitals {x R } centered at the atomic sites {R} 
and specified by the orbital indices {a}. The basis is complete and sufficient to reproduce the one-electron band 
structure of a considered compound in the valent part of the spectrum. The corresponding Hamiltonian matrix is 
denoted by H. The concrete examples of such bases could be the orthonormalized atomic orbitals or the muffin-tin 
orbitals M> We further assume that the basis orbitals are maximally localized (in accordance with a certain criterion 
of maximal localization, the precise form of which is not important here)^ in the sense that any linear combination 
of {xr} will be either less localized or have the same degree of the localization as {xr.}- Then, we immediately 
recognize that the basis functions {xr}) selected in such a way, satisfy to all criteria of WFs, and can be regarded 
as one of possible representations of WFs for the full Hamiltonian H. Indeed, the basis functions {xr} are localized, 
orthonormalized, and (due to the property of completeness of the basis set) any eigenvector of H in the valent part 
of the spectrum can be expressed as a linear combination of {Xr}- This is a natural result and advantage of using 
the localized (or atomic-like) basis. In the plane-wave basis, the localized WFs can be constructed from certain 
number of eigenvectors of H in the valent part of the spectrum, for example, by minimizing the square of the position 
operator (r 2 ). 2 However, we would like emphasize that this is nothing but an elegant way of constructing a compact 
atomic basis from the extended plane waves, a step which becomes rather unnecessary if one works from the very 
beginning in an appropriate basis of atomic orbitals. The precise criterion of the maximal localization (for example, 
the minimization of (r 2 )) is not really important at this stage, because this is merely a mathematical construction 
and depending on the considered physical property one can introduce different criteria of the "maximal localization" . 

However, what we typically need in the process of construction of the model Hamiltonians for the strongly-correlated 
systems is different. For example, for the solution of the many-electron problem it is practically impossible to operate 
with a large number of WFs of the full Hamiltonian H . Instead, one would like to concentrate on the evolution 
of a small number of bands, located near the Fermi level, and construct the WFs only for this group of bands. A 
concrete example is the ti g bands in many transition-metal oxides, which are typically well separated from the rest 
of the spectrum^iS Of course, the WFs for the ti g bands should be orthogonal to the rest of the eigenstates of the 
full Hamiltonian H. This causes an additional complication and the basis functions {Xr}j though can be regarded 
as WFs for the full Hamiltonian H, are no longer the WFs for a limited subspace of eigenstates of H, restricted by 
the tig bands. 

At present, there are two methods, which are commonly used for the construction of WFs for a restricted number 
of states of the full Hamiltonian H: the projector-operator method (for example, employed in Refs. 0, IH, 0, and@) 
and the downfolding method (employed in Refs. IH IToL and ITU) . 
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A. The Projector-Operator Method 



In the projector-operator method, each (nonorthonormalized) Wannier function is generated by projecting a trial 
basis function |x R ) centered at the site R onto a chosen subset of bands (for instance, the t 2g bands, located near the 
Fermi level): 

\W^)=P\xk), (1) 

where 

£=£hMM (2) 

is the projector-operator onto the t^g bands, ipi denotes the eigenstate of H, and i is a joint index combining the 
band index and the position of the momentum k in the first Brillouin zone. The effective tight-binding Hamiltonian 
f)H|f)R.R'll in the region of f 2 g-bands is constructed from the matrix elements of H in the basis of these Wannier 
states: 

C = (Wi\H\W&,). (3) 
The orbitals {W^} shall be further orthonormalized: 

l^) = El<)^ 1/2 ]R'R. ( 4 ) 

R't' 

where <S=||<5>jJ R || is the overlap matrix: 

S& R = {W^W& = {xUP\xk)- (5) 

Then, the effective tight-binding Hamiltonian f) in the basis of orthonormalized Wannier orbitals {1^} takes the 
following form: 

fj = S- 1 ' 2 ^- 1 ' 2 . (6) 

t) is typically regarded as the kinetic-energy part of an effective Hubbard-type model in the projector-operator 
methodi2A£ 

For the subsequent discussions, it will be also convenient to write the overlap matrix S in the form: S—P u , where 
P u is the block of matrix elements of the projector operator @ in the basis of trial orbitals {x R }- Since P commutes 

with H, the projector-operator method guarantees that t) has the same set of eigenvalues in the region of i 2 g-bands 
as the original Hamiltonian H. 



B. The Downfolding Method 



The conventional downfolding method also implies that all atomic basis functions can be divided into two parts: 
{xr}— {Xr}©{Xr}i so that the low-energy part of the spectrum in the direct proximity to the Fermi level is mainly 
composed of the {x R }-states, while {x R } is the rest of the basis states, which mainly contribute to the high-energy 
part. Then, each eigenstate |^) of the full Hamiltonian H can be presented identically as the sum \"4'i} = \'4'i) + \'4>i), 
where \ipf) and \tpl) are expanded over the basis states of correspondingly "i"- and "r"-type, and the matrix equation 
for takes the form: 

(H tt -uM)+H tr \Tpr) = o (7) 
iTU') + ^M) = 0, (8) 

where H 1 ^ 1 ^ are the blocks of matrix elements of H in the basis of "i" ("r")-states. Then, can be expressed 
from Eq. JHJ as 



W)=-(H"- u )- i H rt M) 



(9) 
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and substituted into Eq. J7J. This yields an effective w-dependent Hamiltonian, which acts formally only on \ipl): 

h{cj) = {H u -u)- H tr {H rr - w)- 1 H rt . (10) 

However, \tp\) is only a part of the eigenvector, which is not orthonormalized. Therefore, h(uj) should be additionally 
transformed to an orthonormal basis: 

h(uj) = S- 1/2 {Lu)h(uj)S- 1/2 (cu) + uj. (11) 
This transformation is specified by the overlap matrix 

S{uj) = 1 + H tr (H rr - cu)- 2 H rt , (12) 

which is obtained after the substitution of \ipl), given by Eq. lO, into the normalization condition for the full 

eigenvector \ipi): (ipl\ipl) + (ipl\ipl)—l. h(u>o)t which is typically calculated in the center of gravity of the low-energy 
bands (ljq) constitutes the kinetic-energy part of the effective Hubbard- type model in the downfolding method iriliQr.Un 
Although the downfolding method does not explicitly require the construction of WFs, they are certainly implied 

also in this approach and, at least, can be formally reconstructed from h(u)o)m^ 



C. Downfolding as a Projector-Operator Method 

The conventional downfolding method is exact. However, this property is enforced by the w-dependence of h, 
which is hardly useful from the practical point of view. Formally, for each ipi, oj in Eq. Ijllfl should coincide with 

the eigenvalue of H corresponding to this tpi. Moreover, h{u>) retains an excessive information about H, and the 

full spectrum of eigenvalues and eigenfunctions of the original Hamiltonian H can be formally derived from h(w). 

However, typically we do not need such a redundant information and would like to use h only in order to describe a 
small group of electronic states located near the Fermi level, and do it in the most exact form. 

Therefore, what we want to do next is to show that in a restricted subspace of eigenstates {i/ji} of the original 
Hamiltonian H, the downfolding method can be reformulated as a projector-operator method and retains all attractive 
features of the latter: namely, it becomes exact and does not depend on lu. The trick (which was actually implemented 
m Refs.|3and[l3) is to replace the original Hamiltonian H in the downfolding method by certain modified Hamiltonian 
H, which has the same set of eigenvalues {e^} and eigenfunctions {ipi} in the region of i2g-bands. The rest of the 
eigenstates is not important for our purposes and can be shifted from the valent part of the spectrum to the region 
of infinite energies, specified by the parameter e. Hence, the modified Hamiltonian 7i can be taken in the form: 

u=Yj \i>i) e i{i>i\ + €p ^ = + €p ^ ( 13 ) 

iet 2g 

where P±=l—P is the projector operator to the subspace orthogonal to the i2g-bands. According to the choice of 
the {xr} _ an d {XR.}"basis functions in the downfolding method, the latter mainly contribute to the high-energy part 
of the spectrum, and the overlap between ipi and any linear combination of should be small for the low-energy 

bands. Therefore, all eigenvalues of TrV" r are of the order of e and located in the high-energy part of the spectrum. 
Then, it is intuitively clear that, for e— >oo, the w-dependence in Eq. can be neglected. Therefore, the method 
should be exact and not depend on u>. In order to prove it rigorously, it is convenient to use the idempotency of 
the projector operator: P\=P^. This yields the following identities for the matrix elements p^- r ' t ^ in the basis of 
"t" (V')-states: 

pUptt + ptrprt _ ptt^ ( 14 ) 
pit ptr + ptrprr = ptr ^ ^ 

and 

prtptt + prrprt = prt ( 16 ) 

Then, using Eqs. 114(1. I|15p. and JTBJl, one can prove the following identities, which are valid for the modified Hamil- 
tonian (|13|) in the limit e— >oo: 
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1. The overlap matrix 112[l becomes 

S = (1-P^)- 1 =S~\ 
which is the inverse overlap matrix of the projector-operator method 

2. P" = Pj r (_F[ r ') _1 P£*. Therefore, all terms proportional to e in the Hamiltonian (|lt)|l are exactly cancelled out; 

3. The w-dependent part of the Hamiltonian IjlUI) can be transformed to the form: — ujS^ 1 ; 

4. The remaining part of the Hamiltonian I|1U|) . which depends neither on e nor on oj takes the following form: 

« _ ( J^rr ^ - 1 /jrt „ p*r ^ p.rr ) - 1 ^ rt + ^ ( ) - 1 ^ rr (^rr ^ - 1 ^ ^ where (^t(r)t(r) gtand for the matrix elements of 

the operator PHP in the basis of "f ('V)-states. 

Moreover, it is easy to see that f)" coincides with the Hamiltonian matrix of the projector-operator method prior 
the orthonormalization. Then, by replacing one of the projector operators P in PHP by P 2 and writing explicitly 
the product of two operators in the basis of "t" ("r")-states, one obtains the following identities: 

j^rr prr&rr prt^tr 

jjrt _ prt^tt + prrjjjr* 

and 

£tr _ ^tr prr _|_ j^it ptr 

Using these properties, it is straightforward to verify that the Hamiltonian matrix (|10(l of the downfolding method 
takes the following form: 

h(uj) = S-^S' 1 - US' 1 . 
Then, the orthonormalization transformation yields the following tight-binding Hamiltonian: 

~ h = s-v% tt s- 1 / 2 , 

which is totally equivalent to the tight-binding Hamiltonian (jSJ obtained in the projector-operator method. 

Thus, we have proven that by introducing the modified Hamiltonian of the form l|13|l . the downfolding method can 
be naturally reformulated as the projector-operator method. Therefore, as long as we are interested in the low-energy 
properties of the system, these two methods are equivalent. 

Here, one comment is in order. In our proof, we have used the fact that the overlap between {x R } orbitals and 
the eigenstates {ipi} of the Hamiltonian H should be small for the low-energy bands. However, this implies a specific 
choice for the orbitals {xr} (and also for {x R }), for which the equivalence between the downfolding and the projector- 
operator methods actually takes place. If this property is not satisfied, and the overlap between {x R } and {ipi} is 
large, the downfolding method will eventually breaks down, as we will clearly see it in our calculations for V2O3 
below. 



D. An Analogy with the Many-Electron Problem 

The equivalence between downfolding and projector-operator methods appears to be more generic and has a direct 
implication to the theory of many-electron systems. Indeed, any many-electron state can be expanded over the 
basis of Slater determinants {$}, which play the same role as atomic orbitals in the one-electron version of the 
downfolding method. Suppose that the many-electron ground state is nondegenerate and one can identify a single 
Slater determinant ($g), which has the largest weight in the many-electron wavefunction \^>g) corresponding to the 
ground state. In practice, such \<&g) can be obtained in the frameworks of one-electron Hartree-Fock (or any other) 
theory. Then, each ^ can be formally presented in the form |^ r j) = |^ r *) + |^J'), where is the part proportional 
to |$g) and is expanded over the rest of the basis of Slater determinants. can be further eliminated by 

applying the downfolding method to the true many-electron Hamiltonian H and its eigenfunctions {^i}. By doing 
so, the solution of the many-electron problem can be reformulated in terms of an aj-dependent self-energy E(w)^ in 
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a restricted Hilbert space of many-electron states, which is formally spanned by only one Slater determinant |$g)- 
Nevertheless, such a self-energy contains all information about the ground state as well as the excited states of the 
many-electron system, which is formally incorporated into the w-dependence of S. 

However, the analysis presented in our work clearly shows that if we are interested only in the ground-state properties 
of the many-electron system, the problem can be reformulated in terms of the projection onto the true ground state 
Pg=\^g)(^g\- In this case, the Slater determinant |<£>g) plays the role of the trial wavefunction, which makes 
a formal correspondence between the true many-electron ground state I^g) an d an auxiliary state of independent 
quasiparticles, represented by |$g)- As long as we are concerned with the ground-state properties, such a formulation 
is exact and does not depend on ui. Of course, the direct application of the projector-operator method to the many- 
electron problem is hardly meaningful from the practical point of view, because it implies the knowledge of the true 
many-electron ground state \ ^g) at the first place. However, this simple example, suggested by the projector-operator 
method, clearly shows that, in principle, the procedure of obtaining the total energy of a many-electron system can 
be formulated: 

(a) in a restricted Hilbert space, which is spanned by only one Slater determinant; 

(b) in terms of an w-independent self-energy (or an effective potential) . 

This is a simple but exact property of the total energy of a many-electron system, which does not relies on any 
approximations, ft reminiscences some basic theorems of the density- functional theory (DFT)jiS. and provides some 
insight into why the same many-electron problem can be formulated in two seemingly different ways, one of which is 
DFT and the other one is the self-energy approach. 



III. CHOICE OF TRIAL ORBITALS 

Now let us come back to the construction of WFs and specify the choice of the trial orbitals {x R }- 
As it was already pointed out in the introduction, the basis functions {x R } can be regarded as WFs of the full 
Hamiltonian H. Each basis function is localized around a central atomic site and satisfy certain criterion of the 
"maximal localization" such that any linear combination of {x R } will be "less localized" or at least have the same 
degree of the localization as the original basis set {Xr}- However, this is not necessarily true if one wants to construct 
the WFs only for some part of the electronic structure, specified by certain (restricted) set of eigenstates {4>i} of the 
Hamiltonian H. Due to the additional orthogonality condition to other bands, such a Wannier function will inevitably 
be a linear combination of {xr,} and, hence, a less localized function in comparison with the trial atomic orbital Xr.- 
Nevertheless, one can formulate the problem in a slightly different way and ask which atomic orbital centered at the 
single atomic site is the best representation for the Wannier orbital. Therefore, we search a new set of orthonormalized 
trial orbitals in the form: 

a 

and find the coefficients {c R } from the condition that maximizes the projections (^rJW^c^]), where |VF R [0 R ]) is the 
nonorthonormalized Wannier function constructed from |</> R ) using the projector-operator method: |W / R [<^ R ])=P|0 R ). 
It will automatically guarantee that |0 R ) is the best single-orbital representation for the Wannier function in the 
projector-operator method among the trial orbitals of the form (|1 7|> . By substituting |W R [0 R ]) into expression for 
the projection (0 r |W / r [</>rJ), one can easily find that the problem is reduced to the maximization of the functional 

D = max {(0 R |P|0 R ) - A«&|&> - I)} 

with respect to the coefficient {c R }, where the Lagrange multipliers {A} enforce the orthonormalization condition for 
the new trial orbitals l(T7|) . Then, the maximization of D is equivalent to the diagonalization of Prr=|I(XrI-P|Xr)IIj 
which is nothing but the site-diagonal part of the density matrix constructed for the tig bands in the basis of atomic 
orbitals {xr}- After the diagonalization, we should simply pick up n eigenstates {</> R }, corresponding to the n 
maximal eigenvalues {A}, where n is the number of Wannier orbitals centered at the atomic site Rtf These {</> R } 
will automatically maximize D. This procedure has been proposed in Ref. ||j without proof. Then, some intuitive 
arguments in support of the diagonalization of the density matrix have been given in Ref. lift Here, we have argued 
that it can be derived on the basis of a rigorous variational principle. 
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IV. RESULTS AND DISCUSSIONS 
A. Toy Model 

The simplest model, which nicely illustrates the main idea of the previous Section is the one-dimensional chain 
of atoms {R} (though the spacial dimensionality is not really important here because a similar conclusion can be 
obtained in two and three dimensions). We assume that the atomic basis at each site consists of two orthonormal 
orbitals, say \Xr) an d \Xr)- The full Hamiltonian H is specified by two parameters: the splitting A between atomic 
levels 1 and 2, and the nearest-neighbor transfer integral t (<0), which is assumed to be the same for all pairs of 
orbitals, 

Hrb' = ^ q A ) ^ R ' R ' + ( f t^j {8r,R'+i + $R,R.'-l)- 

The Hamiltonian can be easily diagonalized in the reciprocal (k) space. This yields two eigenstates: \ifiZ) and \tpt)- 
Our goal is to construct the WFs for the lowest eigenstate \ip k ), which can be presented as 

l^fe) =cos9 k \xl) +smd k \x 2 k ), 

where x\ an( i x\ are the Fourier transforms of {x R } an d {Xjj}> respectively, and cos6*fc is nonnegative (that is defined 
by a proper choice of the phase of \tp k ))- If A is large, this eigenstate will be composed mainly by the atomic 
orbitals of the first type, i.e. {xkli which correspond to the largest diagonal matrix elements of the density matrix. 
Therefore, according to the procedure considered in the previous Section, in order to obtain the maximally localized 
representation for the WFs, {xr} should be used as the trial orbitals. 

Now, imagine that we do not know this property and continue to use for the trial orbital \x R ) a linear combination 
of \x R ) and \x 2 R ): 

\Xr)=cosP\xr)+^I3\x r ). 

Then, what will happen? 

It is straightforward to verify that the Wannier function \ W k ~) (in the reciprocal space) obtained after the projection 
of \x R ) onto )(V'fc"l an d the orthonormalization (0J takes the following form 

|MV)=sgn[cos(0 fc -/3)]K7>, 

where sgn[...] stands for the sign of the argument. This means that \WjT) differs from \ipZ) by a phase, which is 
controlled by (3 and depends on k. For small j3, sgn[cos(#fe— /3)]=1 for all k. Therefore, \Wj~) is a smooth function 
in the reciprocal space and its Fourier image is well localized in the real space (Fig.^J. However, when the angle (3 
increases (and starting from certain critical value of (3), sgn[cos(f?fc— (3)] becomes a discontinuous function of k. The 
position of this discontinuity depends both on (3 and on the ratio t/A. Certainly, the discontinuity will affect the 
spread of WFs and make them much less localized in the real space. Thus, the model nicely illustrates the connection 
of the localization of the WFs in the real space with the problem of phase of the Bloch waves in the reciprocal space. 

The degree of localization of the WFs in the real space as a function of (3 is explained in Fig. For /3=0, when 
the trial orbital corresponds to the largest eigenvalue of the local density matrix, the Wannier function spreads only 
over the central and nearest-neighbor atomic sites (which for t/ A=— 0.2 accumulate 99.9% of the total weight of the 
Wannier function). When the angle (3 increases, the weight of the Wannier function transfers (abruptly, starting from 
certain value of (3) from the central part to more remote atomic sites. Hence, the Wannier function becomes less 
localized. This is clearly manifested in the extremely slow convergence of the total weight of the Wannier function in 
the real space (for example, for t/ A=— 0.2 and f3=90°, 99.9% of the total weight can be regained only after including 
about 200 coordination spheres around the central site) . 

Let us discuss the behavior of the model Hamiltonian for the low-energy band in the Wannier basis. The kinetic- 
energy part of the model Hamiltonian JHJl is a quadratic function of WZ ■ Moreover, since H is periodic in the real 

space, its matrix elements will be diagonal with respect to the momentum k in the reciprocal space: \) k =(WjZ\H \WjT). 

Therefore, the phase sgn[cos(0fc— (3)] does not contribute to t) k , and the kinetic-energy part of the model Hamiltonian 
will not depend on the form of the trial orbitals. However, for the Coulomb (and exchange) matrix elements, the 
situation can be different. Although the bare Coulomb interaction is also periodic in the real space, this is a two- 
particle interaction, and its matrix elements in the reciprocal space will contain a combination of the four Wannier 
orbitals: W k _ q W k , +q W k , W k ~ In this case there will be no phase cancelation, and the matrix elements of the Coulomb 
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FIG. 1: (Color online) Results of model calculations for t/A=— 0.2: Wannier functions corresponding to the trial orbitals of 
the form |xk)=cos/3|xfl>+sin/3|x|) with /3=0° (left) and /3=90° (right). The upper panel shows the behavior in the reciprocal 
(k) space, and the lower panel shows the distance-dependence of Wannier functions after the Fourier transformation to the real 
space. The 'Large Component' is the projection onto the atomic orbital \xr), an d the 'Small Component' is the projection 
onto the atomic orbital \x%) (see text for details). 
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FIG. 2: (Color online) Left panel shows the convergence of the total weight of the Wannier function in the real space for two 
types of trial orbitals (specified by the parameter /3). Right panel shows the weight of the Wannier function at the central site 
for two sets of parameters of the model Hamiltonian versus the angle /3. 



interaction will generally depend on the form of trial orbitals. Then, one may ask which representation is better? Of 
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course, all are equivalent and the main question here is whether we would like employ some additional approximations 
or not. A typical approximation is to retain only the site-diagonal part of Coulomb (and exchange) matrix elements 
in the real spaced which may be justified only in the basis of localized WFs (i.e., for /3=0 in the considered example). 
Then, it is reasonable to expect the matrix elements between different lattice sites to be smaller in comparison with 
the site-diagonal part. However, the same approximation is no longer valid for /3=90°, where intersite interactions 
can be of the same order of magnitude as the site-diagonal ones and, therefore, cannot be dropped. 



B. Wannier Functions and Model Hamiltonians for V2O3 

Now let us illustrate how the same procedure works for more realistic systems. For these purposes we pick up rather 
typical (and, in some sense, canonical) example of V2O3 (the space group is D® d , No. 167 in the International Tables). 
The sketch of the crystal structure as well as the positions of the main bands of V2O3 in LDA are summarized in 
Fig. |21 (more detailed information about the lattice parameters and the basis functions used the LMTO calculations 
can be found in Ref . Ho|) . Since each V is surrounded by six oxygen sites, which form a distorted octahedron, the V(3rf) 




FIG. 3: (Color online) Left panel: Total and partial densities of states of V2O3 in the local-density approximation. The shaded 
area shows the contributions of the V(3d)-states. Other symbols show the positions of the main bands. The Fermi level is at 
zero energy. Right panel: the crystal structure of V2O3. 

states are split into the V(3d-t2 g ) and V(3<f-ejp bands, so that the former ones are located near in the low-energy 
part of the spectrum and crossed by the Fermi level within LDA. The trigonal distortion of the VOg octahedra will 
further split the V(3d-t 2g ) levels into one-dimensional a\ g and two-dimensional e g representations. Since the e g and 
e g orbitals belong to the same representation of the point symmetry group, they will mix. The magnitude of this 
mixing is no longer controlled by the symmetry of the system, and can be different for different quantities (depending 
on how these quantities are constructed from the wavefunctions of V2O3). In the following, we will call as "the e g 
orbitals associated with a site R." two eigenstates of the e g symmetry obtained from the diagonalization of the density 
matrix PrrH (XrI-^IXr constructed from the t 2g bands, and corresponding to two largest eigenvalues {A} of Prr- 
Two other eigenstates of the e g symmetry for which {A} are considerably smaller will be called as the e g orbitals (in 
realistic calculations for V2O3, A is about 0.80 and 0.08 for the e g orbitals of the n and a type, respectively). Thus, 
the situation is similar to the toy model considered in Sec. IIV Al 
Then, we take trial orbitals of the e g symmetry in the form 

|XR 9 ) = cos^|0g) + sm/3|g), 

~e" ~e° 

and construct WFs for the t 2g bands. In this expressions, \4>^} and were obtained from the diagonalization 

of the local density matrix, and [3 controls the mixing of the e g orbitals of the 7r and a types. On the contrary, the 
trial orbital of the a\ g symmetry is uniquely determined by the symmetry of the system. Then, according to the 
arguments of Sec. II I II (3=0 should pick up the most localized representation for WFs. This effect is clearly seen in 
our calculations, though there is an important differences from the model considered in Sec. IIV Al Note that in the 
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multi-orbital case, the form of WFs in the reciprocal (k) space is controlled by a matrix [7(k)=||L7 a7 (k)|| (rather than 
a single phase factor) , which generates a new set of WFs after the transformation^ 



\W?) 



cr„ 7 (k)|wz) 



(18) 



In such a case, the change of (3 will generally induce the change of the whole matrix U (k) [for example, through the 
additional orthonormalization transformation ifljl]. which makes some difference from the model analysis presented 
in Sec. IIIII For instance, there will be no cancelation of U(k) in the expression for the matrix elements of H in the 
Wannier basis. Therefore, such matrix elements will generally depend on the representation (3. Nevertheless, we have 
also observed an abrupt change of WFs in the region of large /3, in a close analogy with results of the model analysis 
in Fig. □ 

First, let us discuss the results obtained using the downfolding method for the modified Hamiltonian l|18[l. which 
are shown in Fig. In practical calculations we set e=10 3 Ry. Then, the downfolding method remains stable only 
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FIG. 4: (Color online) The degree of the localization of the Wannier functions for V2O3 depending on the form of trial e g 

orbitals in the downfolding method. The trial orbitals are specified by the transformation )=cos/3|</>p^ )+sin /3|0r ), which 
mixes the e g orbitals of the n and a type. Left panel shows the weights of the Wannier functions at the central V site (the 
Central Weight) and the expectation values of the square of the position operator versus /3. Right panel shows the distance- 
dependence of averaged Hamiltonian f)RR'(d) = (^2 a p ^rr'^r/r) (A being the distance between sites R and R') for two 
values of the parameter j3. 

up to /3=60°. For larger f3, some eigenvalues of Ti rr can fall into the low-energy part of the spectrum, and the whole 
procedure becomes meaningless. The WFs have been reconstructed from the tight-binding Hamiltonian © using the 
method proposed in Ref. flOl which is stable up to /3=45°. 

As the angle (3 increases, the weight of the Wannier function of the e g symmetry at the central V site decreases 
from 0.842 to 0.804 (i.e., by 5%), and the expectation value of r 2 , calculated on these WFs, increases from 2.467 to 
3.183 A 2 (i.e., by nearly 30%). Hence, the WFs become less localized. Is it a big effect or not? In fact, everything 
depends on the considered property. For example, by assuming for a while that all weight of the Wannier function is 
uniformly distributed over the sphere of the radius \J (r 2 ),l& we can estimate (very roughly) the change of the bare 
Coulomb interaction as 12%, which is an appreciable value. 

It is interesting to note that with the increase of /3, the degree of localization of the WFs of the a\ g symmetry also 
decreases, though the effect is considerably smaller. This is not surprising, and the a\ g functions should inevitably 
change because of the additional orthogonalization of them to the e g orbitals. Note that although the Xr 9 and Xr 
orbitals are orthonormal at the same V site, the matrix elements (Xr \P\x'r>) between different atomic sites R and 
R' are generally nonzero. Therefore, the WFs of the a lg and e g symmetry should be additionally orthonormalized, 
and this transformation depends on the value of [3. The same effects leads to the existence of nonvanishing transfer 
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integrals between a\ g and e g orbitals, which are defined as the matrix elements of the tight-binding Hamiltonian f) RR / 
for R^R. 

It is also instructive to consider other possibilities for trial orbitals, which are found from other principles, and 
not directly related with the diagonalization of the density matrix. One possible choice is the e g orbitals corre- 
sponding to the ideal trigonal environment. For the site 1 depicted in Fig. [3] such orbitals have the following formiiS 
(l/v / 3)(|yz)+v / 2|xy)) and (1/V3)(|zx) — V2\x 2 — y 2 )). Another possibility is the so-called "crystal-field orbitals" ob- 
tained from the diagonalization of the site-diagonal part of a more general 5x5 tight-binding Hamiltonian, constructed 
in the basis of all V(3d) orbitals^ Both constructions yield less localized WFs of the e g symmetry in comparison 
with the ones obtained from the diagonalization of the local density matrix. Nevertheless, the difference is very small 
(for example, (r 2 ) = 2.474 and 2.470 A 2 for the ideal trigonal orbitals and the crystal-field orbitals, respectively), 
suggesting that in the case of V2O3, a good starting point for the construction of WFs can be deduced from some 
other (for example, geometrical) consideration. Unfortunately, this rule cannot be applied equally well for other com- 
pounds. For example, one clear exception is the distorted perovskite oxides, where the parameters of the tight-binding 
Hamiltonian in the real space are extremely sensitive to the choice of the local coordinate frame, which is not unique 
and, therefore, the construction of the trial orbitals by means of the diagonalization of the local density matrix seems 
to be the most reliable approachiii 

In the multi-orbital case, the degree of localization of the WFs is related with the distribution of the transfer 
integrals, and it is reasonable to expect that less localized WFs will produce longer-range interactions in the real 
space. This is clearly seen in the behavior of averaged transfer integrals in Fig. If for /3=0° all nonvanishing 
transfer integrals are confined within the radius d<7.5 A(including 16 coordination spheres of the V atoms), those 
for /3=60° have a much longer tail spreading up to d=15 A. At the same time, the energy dispersion obtained after 

the diagonalization of the tight-binding Hamiltonian in the reciprocal space, ik = ER e ~ ,k ' R_R ^RR'i i s practically 
identical for (3= and 60°, and well reproduce the behavior of the original LMTO bands (Fig.[SJ|. 
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FIG. 5: (Color online) LDA energy bands for V2O3 obtained in LMTO calculations and after the tight-binding (TB) 
parametrization using the downfolding method for two values of the parameter (3. Notations of the high-symmetry points 
of the Brillouin zone are taken from Ref. I20L 

Finally, let us discuss the results of the projector-operator method. They are summarized in Fig. El showing the 
weight of WFs at the central V site as a function of (3 and the convergence of the total weight of WFs in the real space 
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for two types of trial orbitals, which are again specified by the parameter (3. When the angle (3 is not particularly 
large so that the downfolding method is applicable, these methods are identical, as it follows from the rigorous proof 
in Sec. Ill Cl For example, the behavior of the central weight in Fig. is practically identical (within the numerical 
accuracy, which is caused by slightly different details of calculations) to that obtained for the downfolding method 
(Fig. |2J|. However, the projector-operator method, which can be supplemented practically with any type of the 
trial orbitals, allows us to consider more extreme scenarios. For example, what may happen if our guess about the 
form of trial orbitals was totally wrong? This situation is nicely illustrated by results of our calculations for /3=90° 
(i.e., by choosing the trial e g orbital in an orthogonal subspace to that, which maximizes local part of the density 
matrix). In this case the Wannier function simply "blows up", and its total weight does not fully converge even 
within d^20 A around the central V siteiSi Thus, the situation is very similar to the model analysis in Sec. II V Al and 
apparently related with the discontinuity of the transformation matrix (|18[) in the reciprocal space. 



1.0 



-a 


0.8* 


Wei 


0.6 




0.4 




0.2 




u 




0.0 







X - a 
• - e 



30 60 
(3 (degrees) 



90 



£ 0.8 X 
"S 0.6 h 




FIG. 6: (Color online) The degree of the localization of the Wannier functions for V2O3 depending on the form of trial e g 

orbitals in the projector-operator method. The trial orbitals are specified by the transformation )=cos/3|0j^ )+sin P\4>-^ ), 
which mixes the e g states of the tt and a type. Left panel shows the weight of the Wannier functions at the central V site (the 
Central Weight) versus (3. Right panel shows the convergence of the total weight of the Wannier functions in the real space for 
/3=90° . The inset shows the amplified data for (3=0°. 



V. SUMMARY AND CONCLUSIONS 



We have considered the problem of construction of the Wannier functions in the first-principles electronic structure 
calculations starting from localized atomic orbitals. There is a number of methods, which are used for these purposes. 
We have argued that some of them can be unified and formulated in an absolutely exact and identical way. We 
have demonstrated this idea for two quite a popular nowadays branches of methods: the so-called projector-operator 
method (considered, for example, in Refs. 0, , and @) and the downfolding method (considered in Refs. 0, ITU 
and ITlf) . We believe that such an analysis will help to resolve a number of controversies existing in this area. Of 
course, identical methods should produce identical results, and we have shown that the root of the differences between 
different applications is related with the choice of the trial orbitals, which are used in order to generate the WFs in 
the projector-operator method or the basis functions for the construction of the effective tight-binding Hamiltonian in 
the downfolding method. We believe that it would be very important in the nearest future to make a similar analysis 
and generalizations for the order- N muffin-tin orbital (NMTO) method*^ However, it is beyond our present abilities. 

The trial orbitals effectively control the localization of WFs. To this end, we have argued that an optimal choice 
for the trial orbitals can be obtained by maximizing the site-diagonal part of the density matrix, constructed from 
the bands of interest. Our numerical calculations show that such a construction already provide a good degree of 
localization for WFs. Of course, depending on the purposes, one can further optimize the WFs (for example, by 
minimizing (r 2 ) or any other quantity). 2,12 However, as a starting point in this procedure one can always use the WFs 
generated by minimizing the site-diagonal part of the density matrix. This is simple and very efficient procedure. We 
have also shown that some deviations from this principle and selection of the trial orbitals in an arbitrary form may 
lead to a substantial derealization of WFs in the real space, which is related with discontinuity of their phase in the 
reciprocal space. 
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